R/Functions to be finished/LeastCost UD/Archive/2019-10-22_lcDistance.R

In vinayudyawer/ATT: The Animal Tracking Toolbox: an extention to the VTrack package

#### Least cost path calculation using passive telemetry datasets
## Written by VU 2018-02-12
## Last updated: 2019-08-29 (VU)
##
## Function to calculate least cost paths between detections from an ATTdata object (from VTrack) 
## around landmasses. Uses 'osmdata' to construct cost raster if none is provided

lcDistance <- function(ATTdata, cost = NULL, utm_epsg, ll_epsg = 4326, timestep = 60, h = 100, UDextent = 1, UDgrid = 100, cost.res = 100, directions = 16, ...){
  ####################################################################################
  ## ATTdata        ATTdata object that consists of Station.Information, Tag.Detections and 
  ##                Tag.Metadata (output from setupData function in VTrack)
  ## cost           cost raster, with land pixel values of 1000 and sea as 1, if none provided
  ##                this is extracted and calculated using polygon downloaded from 'osmdata' package
  ## utm_epsg       EPSG code for CRS object with projected coordinate reference system (units m)
  ## ll_epsg        EPSG code for CRS object with unprojected coordinate reference system (units deg) [default 4326]
  ## timestep       timestep for Center of activity positions (see COA() function in VTrack)
  ## h              smoothing parameter for kernel density estimator in meters (default 100m)
  ## UDextent       a value controlling the extent of the grid used for the estimation of kernel density (default = 4) 
  ## UDgrid         a number giving the size of the grid on which the kernel density should be estimated (in meters; default 100)
  ## cost.res       resolution of cost layer (in degrees) calculated if not provided [default 100]
  ## directions     number of directional axes to consider when calculating 
  ##                least cost path (options: 4, 8, 16 [default])
  ####################################################################################
  
  ## load required libraries and set up CRSs
  sapply(c("lubridate","sf","dplyr","raster","gdistance","spatstat","maptools","rgeos","VTrack", "data.table"), require, character.only=TRUE, warn.conflicts = F, quietly = T)
  ll<-CRS("+init=epsg:4326")
  utm<-CRS(paste0("+init=epsg:", utm_epsg))
  
  combdata <-
    ATTdata$Tag.Detections %>% 
    dplyr::select(-c(Longitude, Latitude)) %>% 
    left_join(ATTdata$Station.Information %>% 
                dplyr::select(Receiver, Station.Name, Station.Latitude, Station.Longitude),
              by=c("Station.Name", "Receiver")) %>% 
    rename(Longitude = Station.Longitude,
           Latitude = Station.Latitude)
  
  statinfo <-
    ATTdata$Station.Information %>% 
    st_as_sf(coords=c("Station.Longitude", "Station.Latitude"), crs=4326)
  
  ## Setup input tagdata; convert to spatial object and transform to projection
  spdata_ll <- 
    combdata %>%
    mutate(step = 0:(nrow(.)-1)) %>% 
    st_as_sf(coords = c("Longitude", "Latitude"), crs = ll_epsg)
  
  spdata_utm <-
    spdata_ll %>% 
    st_transform(crs = utm_epsg)
    
  ## Extract landmass to calculate cost raster if not provided
  if(is.null(cost)){
    message("- No cost layer provided. Downloading coastline data from Open Street Map server")
    require("osmdata")
    tryCatch({
      polydat <-
        opq(bbox = (extent(statinfo) + 0.05)[c(1,3,2,4)]) %>% 
        add_osm_feature(key = 'admin_level', value = '6') %>% 
        osmdata_sf
      
      poly <- 
        polydat$osm_multipolygons %>% 
        st_transform(utm_epsg) %>%
        st_crop(extent(statinfo %>% st_transform(utm_epsg)) + 5000)
      
      cost.ras<-rasterize(poly, raster(extent(statinfo %>% st_transform(utm_epsg)) + 5000, res = cost.res), 1000)
      cost.ras[is.na(values(cost.ras))] <- 1
      projection(cost.ras) <- utm
      
    }, error=function(e){message("Error: no land in sight!\nConsider adding your own cost layer")})
    }else{
      cost.in_utm <- raster::projectRaster(cost, crs=utm, method = "ngb")
      cost.ras <- resample(cost.in_utm, raster(extent(cost.in_utm), res = cost.res), method = "ngb")
      projection(cost.ras) <- utm
      }
  
  ## Produce transition matrices, and correct for distortion
  tryCatch({
    message("- Calculating transition matrices for least cost path estimation")
    trCost <- transition(1/cost.ras, mean, directions = directions)
    trCost <- geoCorrection(trCost, type = "c")
  }, error=function(e){message("Error in calculating Transition layer")})
  
  ## Construct shortest path trajectories between sequence of detection steps
  outdat <- 
    spdata_ll %>% 
    as_Spatial %>% 
    as_tibble %>% 
    rename(Longitude = coords.x1, Latitude = coords.x2) %>% 
    mutate(Distance_m = NA) %>% 
    dplyr::select(-step)
  
  tryCatch({
    message("- Constructing least cost path trajectories between consecutive detections")
    traj<-list()
    for(i in 1:max(spdata_utm$step)){
      if(i %in% 1){pb <- txtProgressBar(min=1, max=max(spdata_utm$step), style=3)}
      
      stepdat <- spdata_utm %>% filter(step %in% c(i-1, i))
      origin <- stepdat %>% slice(1) %>% as_Spatial
      goal <- stepdat %>% slice(n()) %>% as_Spatial
      
      if (nrow(distinct(stepdat)) > 1){
        traj[[i]]<-shortestPath(trCost, origin, goal, output = "SpatialLines")
        outdat$Distance_m[i+1]<-as.numeric(costDistance(trCost, origin, goal))
      } else {
        traj[[i]]<- NULL
        outdat$Distance_m[i+1]<-0
      }
      setTxtProgressBar(pb, i)
    }
    
    traj_clean <- Filter(Negate(is.null), traj)
    trajectory_utm <- do.call(rbind, traj_clean)
    trajectory_ll <- st_as_sf(trajectory_utm) %>% st_transform(ll_epsg)
  }, error=function(e){message("\nError in calculating least cost path trajectories")})
  
  ## Calculate Kernel Density from trajectory
  UDwin <- owin(xrange = (extent(spdata_utm) + (diff(extent(spdata_utm)[1:2])*UDextent))[1:2],
                yrange = (extent(spdata_utm) + (diff(extent(spdata_utm)[3:4])*UDextent))[3:4])
  dimyx <- c(diff(UDwin$yrange)/UDgrid,
             diff(UDwin$xrange)/UDgrid)
  
  # point UD
  tryCatch({ 
    message("\n- Estimating Kernel Density from COA positions")
    coa.data <- 
      COA(ATTdata, timestep = timestep) %>% 
      st_as_sf(coords=c("Longitude.coa","Latitude.coa"), crs=4326) %>% 
      st_transform(utm_epsg) %>% as_Spatial() %>% as_tibble() %>% 
      rename(Longitude = coords.x1, Latitude = coords.x2)
    suppressWarnings(pt.ppp <- ppp(x = coa.data$Longitude, y = coa.data$Latitude, window = UDwin))
    suppressWarnings(UDras.pt <- raster(density(pt.ppp, sigma = h, method = "C", dimyx = dimyx)))
    values(UDras.pt) <- abs(values(UDras.pt)/max(values(UDras.pt)) - 1) * 100
    projection(UDras.pt) <- utm
  }, error=function(e){message("Error in estimating UD from COA positions")})
  
  # trajectory UD
  tryCatch({
    message("- Estimating Kernel Density from least cost path trajectory")
    traj.psp <- as.psp(trajectory_utm, window = UDwin)
    UDras.traj <- raster(density(traj.psp, sigma = h, method = "C", dimyx = dimyx))
    values(UDras.traj) <- abs(values(UDras.traj)/max(values(UDras.traj)) - 1) * 100
    projection(UDras.traj) <- utm
  }, error=function(e){message("Error in estimating UD for trajectories")})
  
  # Merge trajectory UD with point UD and remove land areas
  tryCatch({ 
    UDras <- mean(UDras.traj, UDras.pt)
    mask.ras <- resample(cost.ras, UDras, method = "ngb")
    values(mask.ras)[values(mask.ras) > 1] <- NA
    UDras.m_utm <- raster::mask(UDras, mask = mask.ras)
    UDras.m_ll <- projectRaster(UDras.m_utm, crs = ll)
  }, error=function(e){message("Error in estimating UD")})
    
  # calculate UD areas
  tryCatch({
    message("- Estimating UD area for 50% and 95% contours")
    UD50 <- UDras.m_utm < 50; UD50[values(UD50) %in% 0] <- NA
    UD50.area <- gArea(rasterToPolygons(UD50, dissolve=T))
    
    UD95 <- UDras.m_utm < 95; UD95[values(UD95) %in% 0] <- NA
    UD95.area <- gArea(rasterToPolygons(UD95, dissolve=T))
    
    }, error=function(e){message("Error in estimating UD areas")})
  
  ## return list output with step distances and spatial trajectory file
  out<-list(tagdata = outdat, 
            coa.data = coa.data,
            kernel.areas = data.frame(UD50_m2 = UD50.area, UD95_m2 = UD95.area),
            lc.traj = trajectory_ll, 
            UD.raster = UDras.m_ll, 
            cost.raster = projectRaster(cost.ras, crs = ll, method = "ngb"))
  
  return(out)
}



### Function to plot output on leaflet
plot.lcUD <- function(UDobj, ...){
  library(mapview, warn.conflicts = F, quietly = T)
  library(sf, warn.conflicts = F, quietly = T)
  spdat <-
    UDobj$tagdata %>% 
    st_as_sf(coords=c("Longitude", "Latitude"), crs=4326) %>% 
    group_by(Station.Name) %>% 
    summarise(`Number of Detections` = n())
  
  UDras <- UDobj$UD.raster
  UDras[values(UDras) > 99] <- NA
  
  m <-
    mapview(UDras, layer = "Kernel Utilisation Distribution", homebutton = F, na.color = "transparent", ...) +
    mapview(spdat, alpha = 0, layer = "Detection data", col.regions = "red", cex = "Number of Detections", legend = F, homebutton = F) +
    mapview(st_geometry(UDobj$lc.traj), layer = "Least Cost Path Trajectories", legend = F, homebutton = F)
  
  return(m)
}